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SUMMARY 


We describe a multigrid multiblock method for compressible turbulent flow simulations and present 
results obtained from calculations on a two-element airfoil. A vertex-based spatial discretization 
method and explicit multistage Runge-Kutta time-stepping are used. The slow convergence of a 
single grid method makes the multigrid method, which yields a speed up with a factor of abou 
20, indispensable. The numerical predictions are in good agreement with experimental results. It 
is shown that the convergence of the multigrid process depends considerably on the ordering of the 
various loops. If the block loop is put inside the stage loop the process converges more rapidly than 
if the block loop is situated outside the stage loop in case a three-stage Runge-Kutta. method is usex . 
If a five-stage scheme is used the process does not converge in the latter block ordering. Final y, t e 
process based on the five-stage method is about 60% more efficient than with the three-stage method, 
if the block loop is inside the stage loop. 


INTRODUCTION 


Numerical simulations of turbulent flow in aerodynamic applications are frequently based on th 
Reynolds-averaged Navier-Stokes equations. One of the relevant problems in aeronautics is the pre- 
diction of flow quantities in complicated geometries, such as the multi-element airfoil (see figure 1). 
The simulation of turbulent flow around such a multi-element airfoil configuration was one of the 



Figure 1: Geometry of a two-element airfoil. 

applications selected for the compressible flow solver which was developed by our group and NLR 
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as a part of the Dutch ISNaS project [1]. For this application the use of a single-block, boundary- 
conforming, structured grid is impossible and one may select either an unstructured grid approach 
or a block-structured grid approach. Although the former technique has been successfully applied by 
others [2], we selected the block-structured approach in view of the transparent data structure in the 
coding, ease of implementation of the turbulence model and a high flexibility with respect to the use 
of different physical models in different parts of the computational domain. 

In a previous paper [3] it has been shown that for laminar and turbulent flow around a single airfoil 
the introduction of the multiblock structure has no influence on the results, with respect to both the 
steady-state solution and the convergence rate. Furthermore, invoking the Euler equations instead 
of the Navier-Stokes equations in blocks outside the boundary layer appeared to have no significant 
influence on the results. In this paper we describe the application of the multiblock concept to 
the multi-element airfoil. If the Euler equations are used throughout the computational domain, 
a converged steady-state solution is obtained within a reasonable calculation time. However, if the 
Reynolds-averaged Navier-Stokes equations are solved in the boundary layers, the rate of convergence 
is unacceptably low. Therefore, a multigrid technique was implemented in order to accelerate the 
convergence. The resulting gain in calculation time is close to a factor of 20, and the converged 
solution is in good agreement with wind-tunnel measurements. 

In section 2 the numerical technique, which is based on a combination of a finite volume method 
with central spatial differencing and a Runge-Kutta explicit time-stepping method, is described. The 
results, both for inviscid and for viscous simulations, are presented in section 3. Finally, in section 4 
some conclusions are summarized. 


NUMERICAL METHOD 


In this section we describe the numerical method used in the flow solver. The two-dimensional, 
compressible Navier-Stokes equations can be written in integral form as 


d_ 

dt 




(i) 


where U represents the vector of dependent variables, 


U = [p, pu, pv, E] t , (2) 

with p the density, u and v the Cartesian velocity components, and E the total energy density. 
Further, O is an arbitrary part of the two-dimensional space with boundary OQ and F and G are 
the Cartesian components of the total flux vector. This flux vector consists of two parts: the non- 
dissipative or ’convective’ part and the dissipative or ’viscous’ part, which describes the effects of vis- 
cosity and heat conduction, and involves first order spatial derivatives. The Navier-Stokes equations 
(1) are averaged over a sufficiently large time interval. Due to the nonlinear terms in the convective 
fluxes, the resulting ’Reynolds-averaged Navier-Stokes’ equations involve averages of products of two 
velocity components. These terms are modeled by a suitable turbulence model. In the present paper 
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the algebraic Baldwin-Lomax turbulence model, in which the unknown terms are modeled by eddy 
viscosity terms, is adopted [4]. 

The discretization of the Navier-Stokes equations follows the method of lines, i.e. the spatial 
discretization is performed first, and subsequently the resulting set of ordinary differential equations 
is integrated in time, until the steady state solution is approximated. First the computational domain 
is divided into blocks and each block is partitioned in quadrilateral cells with the help of a structured, 
boundary-conforming grid. The variables are stored in the grid points. A finite volume method is used 
in which the integral form of the Navier-Stokes equations is applied to a control volume ft, bounded 
by the dashed lines in figure 2. The convective flux through a boundary of this control volume is 



Figure 2: Control volume in the vertex-based method. 

approximated using the value of the convective flux vector in the midpoint of the boundary. The 
latter is calculated by averaging over the two neighboring grid points. The viscous flux vector involves 
spatial derivatives of the state vector U and is approximated in the corner points of the control volume 
with the use of Gauss’ theorem on a grid cell. The viscous flux is subsequently calculated using the 
trapezoidal rule. This method is called the vertex-based method. 

The method of central differencing leads to a decoupling of odd and even grid points and to 
oscillations near shock waves. Even in viscous flow calculations the presence of the viscous dissipation 
is insufficient to damp these instabilities outside shear layers. Therefore, nonlinear artificial dissipation 
is added to the basic numerical scheme. This artificial dissipation consists of two contributions: fourth 
order difference terms which prevent odd-even decoupling, and second order difference terms to resolve 
shock waves. The second order terms are controlled by a shock sensor, which detects discontinuities 
in the pressure. In the present flow solver the artificial dissipation in the boundary layers, where 
the viscous dissipation should be dominant, may be reduced by multiplication with the ratio of the 
local and free-stream Mach number. The role of the artificial dissipation in relation to the viscous 
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dissipation is discussed in more detail in reference [5], 

At the solid wall boundaries the no-slip condition is used. The density and energy density in the 
grid points on a solid wall are calculated by solving the corresponding discrete conservation laws, 
using the two adjacent cells within the computational domain and their mirror images inside the wall 
as the control volume. The values of the density and energy density in the grid points inside .the walls 
are adjusted such that the adiabatic wall condition is approximated. The boundary conditions at 
a (subsonic) far- field boundary are based on characteristic theory. The extent of the computational 
domain can be reduced without affecting the accuracy if a vortex is superimposed on the incoming 
free stream outside the computational domain [6]. 

Due to the topology of the two-element airfoil geometry, special points in the computational grid are 
unavoidable. The computational grids used contain two special points at block boundaries, where five 
cells meet (see figure 4). These points can be treated in an elegant way within the same numerical 
scheme, if the dummy vertices outside the ’current 7 block are defined appropriately. The multi- 
valuedness of the variables at the special point, caused by this asymmetric treatment, is eliminated 
by taking the average of the five different values after all blocks have been treated. This is sketched 
in figure 3. 



Figure 3: Control volume for a special point. 


The system of ordinary differential equations, which results after spatial discretization, is integrated 
in time using a time-explicit multistage Runge-Kutta method. In the present flow solver a three-stage 
scheme in which the dissipative fluxes (both viscous and artificial) are calculated once per time-step, 
and a five-stage scheme in which the dissipative terms are calculated only at the odd stages, are 
implemented. With this treatment both calculation time is saved and the stability region of the 
method is increased. Extra calculation time is saved by advancing each grid point at the maximum 
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local time-step according to its own stability limit. In this way the evolution from the initial solution 
to the steady state is no longer time accurate, but the steady state solution obtained is unaffected. 

The above time-stepping method acts as the relaxation method and coarse grid operator in the 
multigrid solver (see reference [6]). In this solver an initial solution on the finest grid is obtained with 
a full multigrid method. This initial solution is corrected in the FAS-stage, where either V- or W- 
cycles can be chosen. A fixed number of pre- and post-relaxations is performed before turning to the 
next coarser or finer grid. The solution is transferred to a coarser grid by injection, the residuals by 
full weighting and the corrections to the solution are prolonged by bilinear interpolation. In order to 
increase the smoothing properties of the Runge-Kutta time-stepping technique an implicit averaging 
of the residuals is applied with frozen residuals at the block boundaries. For mono-block applications 
this method has given satisfactory results for both two-dimensional and three-dimensional flows [5]. 

In the multi-element airfoil application care has to be taken in the definition of the residual- vector 
in the special points. The proposed treatment of a special point implies that the control volume is 
different in each of the five blocks where such a point is found. In the required averaging the five 
residual-vectors in a special point are weighed with their corresponding time-steps. Without this 
weighing the multigrid process cannot converge to the single grid stationary state solution. 

In this multigrid, multiblock solver with a multistage time-stepping method there are various 
possibilities for intertwining the different loops. In the present study the grid loop is chosen as the 
outer loop and the effect of interchanging the block and the stage loop will be studied. Several 
’competing’ requirements serve as possible guidance for selecting a specific ordering of these loops. 
On the one hand an anticipated parallel processing of the different blocks is more efficient, if the 
data transfer between the blocks is kept to a minimum, i.e. with the stage loop inside the block 
loop. On the other hand the good convergence of the multigrid mono-block solver may be reduced as 
the dummy variables near the block boundaries are kept frozen during more stages of the time-step. 
This would suggest to put the block loop inside the stage loop. In order to study this dilemma we 
implemented these two loop orders in a flexible way: a single parameter determines whether the block 
loop is situated inside or outside the stage loop. 


RESULTS 


Description of the test-case 


We will present results for a two-component airfoil geometry consistii^g of the NLR7301 wing 
section, from which a flap has been cut out at a deflection angle of 20° and with a gap width of 2.6% 
chord length [7] (see figure 1). The combination of a Mach number of 0.185 and an angle of incidence 
of 6° or 13.1°, of which the latter is close to maximum lift conditions, yields subsonic flow. The 
Reynolds number based on the chord length of the airfoil is 2.51 x 10 6 . In the viscous calculations 
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the locations of the transition from laminar to turbulent flow are prescribed. 

The C-type computational grids (either for inviscid or viscous flow) were constructed by J. J. Benton 
from British Aerospace, and are subdivided in 37 blocks (see figure 4). The grid lines are continuous 
over block boundaries. Two grids are used: one ’Euler’ grid (inviscid) consisting of 16448 cells, and 
a ’Navier-Stokes’ grid (viscous), which is refined in the boundary layers and wakes and consists of 
28288 cells. 



Figure 4: Block structure of the computational grid. 


For both angles of incidence results from wind-tunnel measurement by Van den Berg [7] are avail- 
able, including velocity profiles in the boundary layers and the pressure coefficient on the profile. 
Since the flow is attached apart from a small laminar separation bubble near the leading edge of 
the wing, the adopted turbulence model should be adequate and yield a useful comparison between 
experiment and calculation. 


Inviscid Flow 


In order to test the flow solver on the complicated block structure of the two-element airfoil geom- 
etry, we considered the relatively simple inviscid flow case, where in all blocks the Euler equations are 
solved. In this way problems related to the turbulence model are separated from possible algorithmic 
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problems. The use of the Euler equations implies that the boundary conditions at the solid wall 
boundaries have to be changed. For inviscid flow there is only one physical boundary condition of 
zero mass flux through the wall. In the vertex based approach the density, the pressure and the 
tangential velocity at the wall are approximated by linear extrapolation. 

In figure 5 the multigrid convergence behavior of the solver in the 13.1° case is shown. The discrete 
L 2 -norm of the residual of the density is plotted as a function of the number of W-cycles. A converged 
solution is obtained within a much smaller calculation time when compared to the single grid approach 
even though only three different grid levels are available. Both for the single grid and the multigrid 
calculations machine accuracy was obtained. The specific block structure nor the treatment of the 
special points leads to any specific difficulties. For this inviscid test a comparison with experimental 
results is not meaningful and will not be made. 



Figure 5: Convergence behavior for inviscid flow at an angle of incidence of 1 3 . 1 c 


Viscous Flow 


We consider the simulations of turbulent, viscous flow and present results for the 6° case only. 
Si'hgle-grid calculations in which only local time-stepping is applied as a convergence acceleration 
technique yield a steady-state solution which is in good agreement with the experimental results. 
However, in contrast with a fully inviscid simulation, the rate of convergence is very small, and 
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renders this method unacceptable for practical applications. Therefore, as a method to increase 
the convergence rate further, the multigrid technique and implicit residual averaging as described in 
section 2 are indispensable. 

In a simulation of turbulent flow at high Reynolds number it is important that the effects related 
to the physical dissipation are not outweighed by those of the numerical or artificial dissipation. This 
requirement could give rise to difficulties in the present multigrid method, since the time-stepping 
method used requires a certain minimum amount of dissipation for sufficient smoothing of the large 
wave-number components of the error (see reference [5]). If the artificial dissipation in the boundary 
layer is reduced by scaling with the ratio of the local and free-stream Mach number, i.e. decreasing 
the smoothing properties of the time-stepping method, a converged solution (engineering accuracy) 
could be obtained by increasing the number of pre- and post-relaxations. The convergence behavior of 
this calculation during the FAS stage is shown in figure 6, where the discrete X 2 -norm of the residual 
of the density is plotted as a function of the number of W-cycles. In the blocks outside the boundary 
layers and wakes the Euler equations are solved instead of the Navier-Stokes equations. The good 



Figure 6: Viscous flow at an angle of incidence of 6.0°: convergence behavior 

agreement with the wind-tunnel measurements can be inferred from figure 7, where the experimental 
and numerically predicted pressure coefficients on the airfoil and flap are shown. 

This solution was obtained with the block loop inside the stage loop of the five-stage Runge-Kutta 
time-stepping method. Hence, the variables at the dummy vertices outside a block are updated 
after every stage, which implies that the effects of the multiblock structure on the convergence are 
kept to a minimum. The frequency of data transfer between the blocks makes this method less 
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Figure 7: Viscous flow at an angle of incidence of 6.0°: comparison of the pressure coefficient on the 
airfoil between calculation (solid) and experiment (dashed). 

efficient for parallel processing. However, with the block loop outside the stage loop, i.e. with an 
update of the dummy variables only after five flux evaluations, a converged solution could not be 
obtained. Apparently, the interval between two moments of data transfer between the blocks has to 
be sufficiently small in order to obtain a convergent multigrid method. 

Further evidence for this statement is obtained from calculations with a three-stage instead of a 
five-stage Runge-Kutta time-stepping method. If the block loop is outside the stage loop, the dummy 
variables are updated after three flux evaluations. Although the rate of convergence is lower than 
in the case with the loops interchanged (see figure 8), the solution has converged within engineering 
accuracy after ~ 200 W-cycles. A comparison of the three-stage and five-stage schemes with the 
block loop inside the stage loop shows that the five-stage scheme is more efficient: about 60 W-cycles 
suffice to get the residuals at the same level as with the three-stage scheme after 200 W-cycles. The 
five-stage scheme leads to a reduction in calculation time of approximately 60% in this instance. 


DISCUSSION 


We presented simulation results obtained with a multigrid multiblock method for a two-element 
airfoil. Both viscous and inviscid calculations were performed using the same multigrid process 
and the same vertex-based spatial discretization method. Moreover, either a three- or a five-stage 
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Figure 8: Convergence behavior of the three-stage Runge-Kutta scheme for turbulent flow; comparison 
between block loop inside (solid) and outside (dashed) stage loop. 

Runge-Kutta scheme was considered for the integration in time and the smoothing properties of this 
relaxation method were further enhanced through the introduction of local time-stepping, implicit 
lesidual averaging in which the residuals at the block boundaries were kept fixed to their non-smoothed 
values. 

The mviscid calculations have shown that a solution which is converged up to machine accuracy can 
be obtained with this multigrid method. A comparison with the single grid simulation method shows 
that a considerable reduction in calculation time was obtained with the multigrid method, although 
the convergence of the single grid method for inviscid calculations was already quite acceptable. We 
also investigated two different numerical boundary conditions at the solid walls. It appeared that 
linear extrapolation of the pressure not only leads to a better convergence than constant extrapolation, 
but also gives rise to a much smaller entropy layer around the airfoil. The resulting drag coefficient, 
which theoretically should equal zero in this subsonic flow, is reduced by almost 60%. 

In the viscous calculations the single grid method was found to yield a well converged result in the 
6 -case, however, the convergence towards the steady state solution was extremely slow and makes 
the use of a multigrid approach essential. A comparison of the calculation times required in both 
methods shows that a total speed-up with a factor of about 20 can be reached. The numerical 
predictions obtained for the lift- and pressure coefficients compare well with experimental results 
and give confidence in the use of the Baldwin-Lomax model for this application. The convergence 
of the multigiid process was studied in detail, showing that the ordering of the various loops in the 
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process has a considerable effect. Interchanging the block and stage loops and keeping the grid loop as 
the outer loop, yields an optimal convergence when the block loop is put inside the stage loop. If the 
stage loop is put inside the block loop then convergence of the multigrid process was absent when 
using the five-stage Runge-Kutta method as the relaxation method. Apparently, the smoothing of 
the relaxation method becomes less effective as the number of stages between two ’updates’ of the 
dummy-variables increases. This result has some less favorable consequences in view of a possible 
parallel processing of the multigrid method. On the one hand parallel processing seems more efficient 
if the frequency of data transfer between the blocks can be reduced. On the other hand the reduction 
of this frequency results in a reduction of the convergence rate of the multigrid process, and in some 
instances even to an absence of convergence. This suggests that in a possible parallel processing of 
this multigrid method, an optimal rate of data-exchange between the blocks should be determined. 
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